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£2 : ABSTRACT 

We present results of a convergence study in which we compare the density profiles of CDM 
dark matter halos simulated with varying mass and force resolution. We show that although 
increasing the mass and force resolution allows one to probe deeper into the inner halo regions, 
the halo profiles converge at scales larger than the "effective" spatial resolution of the simulation. 
This resolution is defined both by the force softening and by the mass resolution. On radii 
larger than the "effective" spatial resolution, density profiles do not experience any systematical 
trends when the number of particles or the force resolution increase further. In the simulations 
presented in this paper, we are able to probe density profile of a relaxed isolated galaxy-size halo 
at scales r = (0.005 — l)r v i r . We find that the density distribution at resolved scales can be 
well approximated by the profile suggested by Moore et al.(I998): p oc x~ 15 (l + x 1,5 ) -1 , where 



x = r/r s and r s is the characteristic radius. The analytical profile proposed by Navarro et al. 
(1996) p oc a: _1 (l + x)~ 2 , also provides a good fit, with the same relative errors of about 10% for 
radii larger than 1% of the virial radius. For this limit both analytical profiles fit well because 
for high-concentration galaxy-size halos the differences between these profiles become significant 
only at scales well below O.Olrvir- We also find that halos of similar mass may have somewhat 
different parameters (characteristic radius, maximum rotation velocity, etc.) and shapes of their 
density profiles. We associate this scatter in properties with differences in halo merger histories 
and the amount of substructure present in the analyzed halos. 

Subject headings: cosmology:theory - galaxy structure - methods: numerical 

1. Introduction nated dwarf galaxies are at odds with predictions 

of hierarchical models. Specifically, it was argued 
During the last decade there has been an in- (Flores & Primack 1994; Moore 1994) that circu- 

creasingly growing interest in testing the predic- lar ve l cities, v c (r) = [GM (< r)^] 1 ' 2 , at small 

tions of variants of cold dark matter (CDM) mod- galactocentric radii predicted by the models are 

els at subgalactic (< 100 kpc) scales. This interest too m g n an( j mcrease too rapidly with increasing 

was initiated by indications that observed rotation radius compar ed to the observed rotation curves, 

curves in the central regions of dark matter domi- The ste eper than expected rise of v c {r) implies 

" that the shape of the predicted halo density dis- 

tribution is incorrect and/or that the DM halos 
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formed in CDM models are too concentrated (i.e., 
have too much of their mass concentrated in the 
inner regions). 

In addition to the density profiles, there is an 
alarming mismatch in the predicted abundance of 
small-mass (< 10 8 — 10 9 /i~ 1 M Q ) galactic satellites 
and the observed number of satellites in the Lo- 
cal Group (Kauffmann, White & Guiderdoni 1993; 
Klypin et al. 1999; Moore et al. 1999). Although 
this discrepancy may well be due to feedback pro- 
cesses (such as photoionization) which prevent gas 
collapse and star formation in the majority of the 
small- mass satellites (e.g., Bullock, Kravtsov & 
Weinberg 2000), the mass scale at which the prob- 
lem sets in is similar to the scale in the spectrum 
of primordial fluctuations that may be respon- 
sible for the problems with density profiles. In 
the age of precision cosmology that forthcoming 
MAP and Planck cosmic microwave background 
anisotropy satellite missions are expected to bring, 
tests of the cosmological models at small scales 
may prove to be the final frontier and the ulti- 
mate challenge to our understanding of cosmology 
and structure formation in the Universe. How- 
ever, this obviously requires detailed predictions 
and checks from the theoretical side, as well as 
higher resolution/quality observations and a good 
understanding of their implications and associated 
caveats. In this paper we focus on the theoretical 
predictions of the density distribution of DM ha- 
los. 

A systematic study of halo density profiles for 
a wide range of halo masses and cosmologies was 
done by Navarro, Frenk & White (1996, 1997; 
hereafter NFW), who argued that the analytical 
profile of the form p(r) = p s (r/r s ) _1 (l +r/r s )~ 2 
provided a good description of halo profiles in 
their simulations for all halo masses and in all 
cosmologies. Here, r s is the scale radius which, 
for this profile, corresponds to the scale at which 
d log p(r)/d log r | r=rs = —2. The parameters of 
the profile are determined by the halo's virial mass 
M v i r and concentration defined as c e r w - 1T /r s . 
NFW argued that there is a tight correlation be- 
tween c and M v ; r , which implies that the density 
distributions of halos of different masses can in fact 
be described by a one-parameter family of analyt- 
ical profiles. Further studies by Kravtsov, Klypin 
& Khokhlov (1997), Kravtsov et al. (1998, here- 
after KKBP98), Jing (2000), Bullock et al. (2000), 



although confirming the c(M v ; r ) correlation, indi- 
cated that there is significant scatter in both the 
density profiles and concentrations for DM halos 
of a given mass. 

Following the initial studies by Flores & Pri- 
mack (1994) and Moore (1994), KKBP98 pre- 
sented a systematic comparison of the results of 
numerical simulations with rotation curves of a 
sample of seventeen dark matter dominated dwarf 
and low surface brightness (LSB) galaxies. 

We pointed out that the measured rotation 
curves of these galaxies all had the same shape 
with nearly linear central behavior, and further- 
more, based on comparison with the density pro- 
files of simulated halos there did not seem to be 
a significant discrepancy in the shape of the den- 
sity profiles at the scales probed by the numer- 
ical simulations (> 0.02 — 0.03r v i r , where r v - lr is 
halo's virial radius). In other words, the central 
density distribution in both galaxies and CDM ha- 
los was found to be shallower than oc r _1 . These 
conclusions were subject to several caveats and re- 
quired further testing. First, observed galactic ro- 
tation curves had to be re-examined more care- 
fully and with higher resolution. The fact that 
all of the observed rotation curves used in earlier 
analyses were obtained using relatively low resolu- 
tion HI observations required checks of the possi- 
ble beam smearing effects. Also, the possibility of 
non-circular random motions in the central regions 
which could modify the rotation velocity of the 
gas (e.g., Binney & Tremain 1987, p. 198) had to 
be considered. Second, the theoretical predictions 
had to be tested for convergence and extended to 
scales < 0.01r v i r . 

Moore et al. (1998; see also a more recent 
convergence study by Ghigna et al. 1999) pre- 
sented a convergence study arguing that mass res- 
olution has a significant impact on the central den- 
sity distribution of halos. They suggested that 
at least several million particles per halo are re- 
quired to reliably model the density profiles at 
scales < 0.01r v ; r . Based on these results, Moore 
et al. (1998) advocated a density profile of the 
form p(r) oc (r/r s ) _1 - 5 [l + (r/rg) 1 - 5 ] -1 , that be- 
haves similarly (p oc r -3 ) to the NFW profile at 
large radii, but is steeper at small r: p a r 1,5 . 
Most recently, Jing & Suto (2000) presented a sys- 
tematic study of density profiles for halo masses in 
the range 2 x IO^/i^Mq - 5 x 10 14 ft- 1 M Q . The 
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study was uniform in mass and force resolution 
featuring ~ 5 — 10 x 10 5 particles per halo and 
force resolution of w 0.004r V j r . They found that 
galaxy-mass halos in their simulations are well fit- 
ted by the profile 2 p(r) oc (r/r s )~ 1 - 5 [l + r/r s ]~ 1 - 5 , 
but that cluster-mass halos are well described by 
the NFW profile, with logarithmic slope of the 
density profiles at r = 0.01r v ; r changing from 
w -1.5 for M vir ~ 10 12 /i _1 M Q to « -1.1 for 
M vir - 5 x 1O 14 /i _1 M . Jing & Suto interpreted 
these results as an evidence that profiles of DM 
halos are not universal (but see §3.1 for a possible 
alternative interpretation). 

At small scales, the results of Kravtsov et al. 
(1998) are at odds with the results of above stud- 
ies. Although fairly extensive convergence tests 
were done in that study, they focused on the ef- 
fects of spatial resolution and the mass resolution 
was kept constant in almost all the tests. In this 
case we found that halo density profiles converged 
at scales larger than two formal resolutions of the 
ART code. As we will show in this paper, this 
is not true when the mass resolution is varied. In 
particular, the convergence study presented in this 
paper in which we varied both mass and force reso- 
lution shows that for the ART simulations conver- 
gence is reached at scales larger than four formal 
resolutions or the scales containing 200 particles, 
whichever is larger. The shallow behavior of the 
density profiles in Kravtsov et al. was found at 
scales w 2 — 5 formal resolutions (at larger scales 
profiles were consistent with the NFW functional 
form) and is therefore a numerical artifact. In sim- 
ulations presented in this paper we find profiles 
that are consistent with cuspy NFW and Moore 
et al. distributions at well resolved scales. 

New observational and theoretical develop- 
ments show that comparison between model pre- 
dictions and observational data is not straight- 
forward. Decisive comparisons require reaching 
convergence of theoretical predictions and under- 
standing the kinematics of the gas in the central 
regions of observed galaxies. As we noted above, in 
this paper we present convergence tests designed 
to test effects of mass resolution on the density 

2 Note that this profile is somewhat different than the profile 
advocated by Moore et al., but behaves similarly to the 
latter at small radii. Figure 9 shows that all three profiles 
— NFW, Moore, and Jing & Suto — provide good fits to 
dark matter halos simulated at high resolution. 



profiles of halos formed in the currently popular 
CDM model with cosmological constant (ACDM) 
and simulated using the multiple mass resolu- 
tion version of the Adaptive Refinement Tree code 
(ART). This study is crucial in resolving the dis- 
crepancy of our previous study on the density 
profiles with other numerical studies. We also 
discuss several caveats with respect to drawing 
conclusions about the density profiles from the 
fits of analytical functions to numerical results 
and their comparisons to observational data. In 
the following section we describe the code and 
numerical simulations used in our analysis. In §3 
we compare the analytical fits advocated by NFW 
and Moore et al., fits of these profiles to the den- 
sity profiles of simulated halos, and convergence 
analysis of our numerical results. 

2. Numerical simulations 
2.1. Code description 

The Adaptive Refinement Tree code (ART; 
Kravtsov, Klypin & Khokhlov 1997) was used to 
run the simulations. The ART code starts with 
a uniform grid, which covers the whole computa- 
tional box. This grid defines the lowest (zeroth) 
level of resolution of the simulation. The stan- 
dard Particles-Mesh algorithms are used to com- 
pute density and gravitational potential on the 
zeroth-level mesh. The ART code reaches high 
force resolution by refining all high density regions 
using an automated refinement algorithm. The 
refinements are recursive: the refined regions can 
also be refined, each subsequent refinement hav- 
ing half of the previous level's cell size. This cre- 
ates a hierarchy of refinement meshes of different 
resolution, size, and geometry covering regions of 
interest. Because each individual cubic cell can 
be refined, the shape of the refinement mesh can 
be arbitrary and match effectively the geometry of 
the region of interest. 

The criterion for refinement is the local density 
of particles: if the number of particles in a mesh 
cell (as estimated by the Cloud-In-Cell method) 
exceeds the level n t hresh, the cell is split ( "refined" ) 
into 8 cells of the next refinement level. The re- 
finement threshold may depend on the refinement 
level. The code uses the expansion parameter a 
as the time variable. During the integration, spa- 
tial refinement is accompanied by temporal refine- 
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ment. Namely, each level of refinement, I, is in- 
tegrated with its own time step Aai = Aao/2 l , 
where Aao is the global time step of the zeroth 
refinement level. This variable time stepping is 
very important for accuracy of the results. As the 
force resolution increases, more steps are needed 
to integrate the trajectories accurately. In the re- 
mainder of the paper by the term formal resolu- 
tion, /iformai, we will mean the size of a cell on the 
highest level of refinement reached in simulation. 
This is similar to the usual practice of defining 
formal resolution in the uniform grid codes. The 
actual force resolution of the code is somewhat 
larger than the formal resolution. In the ART 
code, the interparticle force is weaker ("softer") 
than the Newtonian force at scales < 2/if orma i. The 
average interparticle force is Newtonian at scales 
^ 2/iformai; although there is a substantial scatter 
in the force at w 2 — 3/if or mai due primarily to er- 
rors in numerical differentiation of potential. For 
comparison, average Plummer softened force with 
softening epiummcr = ^formal, reaches the Newto- 
nian value at scales « 5 — 6h{ oima ,\ (see Kravtsov et 
al. 1997 for details). Therefore, one formal reso- 
lution in the ART code is equivalent to « 0.3 — 0.4 
Plummer softening of the same value. 

The cosmological simulations performed using 
the ART code were compared with the simula- 
tions started from identical initial conditions and 
performed using the well-known PM and AP 3 M 
codes. The comparisons showed that results of 
ART simulations (for a wide battery of the com- 
monly used statistics and halo parameters) are 
similar to those of the AP 3 M simulations at all re- 
solved scales. These comparisons and other tests 
of the ART code can be found in Kravtsov (1999) 
and Knebe et al. (2000). 

2.2. Initial conditions 

The current version of the ART code has the 
ability to handle particles of different masses. In 
the present analysis this ability was used to in- 
crease the mass (and correspondingly the force) 
resolution inside a few pre-selected halos. The 
multiple mass resolution is implemented in the fol- 
lowing way. We set up a realization of the ini- 
tial spectrum of perturbations in such a way that 
a very large number of small-mass particles can 
be generated in the simulation box. For example, 
for the first (second) set of simulations (see be- 



low) we generate 512 3 (1024 3 ) independent spec- 
trum harmonics. Potentially, initial conditions 
with 512 3 /1024 3 particles could be generated. Co- 
ordinates and velocities of the particles are cal- 
culated using all waves ranging from the funda- 
mental mode k = 2-k/L to the Nyquist frequency 
k = 2n/L x A fl / 3 /2, where L is the box size and 
TV is the number of particles in the simulation. 

The code actually generates positions and ve- 
locities for all 512 3 /1024 3 particles, but some of 
the particles are then merged into particles of 
larger mass. The larger mass (merged) particle 
is assigned velocity and displacement equal to the 
average velocity and displacement of the merged 
particles. The whole lagrangian space of particles 
is divided into large cubic blocks of particles with 
each block having 16 3 particles. Depending on 
what local mass resolution is required, each par- 
ticular block can be subdivided into smaller sub- 
blocks and generate from 1 to 16 3 particles (the 
highest resolution). Using this procedure 3 , we can 
generate particles with 5 different masses covering 
dynamic mass range of 4096. 

We start simulations by making a low resolu- 
tion run with uniform mass resolution in which all 
particles have the largest possible mass. Next we 
run simulations with 32 3 and 64 3 particles. Using 
these runs, we identify halos in the simulation and 
select halos to be re-simulated with higher mass 
and force resolution. For each selected halo we de- 
termine its virial radius. We then identify all par- 
ticles inside the two virial radii and find lagrangian 
coordinates of each particle. The coordinates are 
used to mark blocks of particles to generate the 
initial conditions of the highest mass resolution. 
Once all particles are processed and all blocks are 
marked, we mark all blocks adjacent to those al- 
ready marked to produce initial conditions of the 
eight times lower mass resolution. This procedure 
is repeated for lower and lower mass resolution lev- 
els. In the end, each unmarked block will produce 
one most massive particle and a marked block will 
generate a number of particles which depends on 
the step in which the block was marked. Figure 1 
shows the outcome of the process of mass refine- 
ment in a 2-dimcnsional case. 



3 The code is actually written to handle an arbitrary dynamic 
range. The current limit is determined by computational 
limitations. 
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Fig. 1. — Example of the construction of mass re- 
finement in lagrangian space (here for illustration we 
show a 2D case). Three central blocks of particles 
were marked for highest mass resolution. Each block 
produces 16 2 particles of the smallest mass. Adjacent 
blocks correspond to the four times lower resolution 
and produce 8 2 particles each. The procedure is re- 
peated recursively until we reach the lowest level of 
resolution. The region of the highest resolution can 
have arbitrary shape. 

Figure 2 shows an example of mass refinement 
for one of the halos in our simulations. A large 
fraction of high resolution particles ends up in 
the central halo, which docs not have any larger 
mass particles (see insert in the bottom panel). At 
z = 10, the region occupied by the high resolution 
particles is non-spherical: it is substantially elon- 
gated in the direction perpendicular to the large 
filament clearly seen at z = 0. 

After the initial conditions are set, we run the 
simulation again allowing the code to perform 
mesh refinement based only on the number of par- 
ticles with the smallest mass. 

2.3. Numerical simulations 

We simulated a flat low-density cosmological 
model (ACDM) with Q a = 1 - Q A = 0.3, the 
Hubble parameter (in units of 100 kms _1 Mpc _1 ) 
h = 0.7, and the spectrum normalization erg = 0.9. 
We have run two sets of simulations. The first set 
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Fig. 2. — Distribution of particles of different masses 
in a thin slice through the center of halo Ai (see Ta- 
ble 1) at z = 10 (top panel) and at z = (bot- 
tom panel). To avoid crowding of points the thick- 
ness of the slice is made smaller in the center (about 
30/i _1 kpc) and larger (lfe _1 Mpc) in the outer parts 
of the forming halo. Particles of different mass are 
shown with different symbols: tiny dots, dots, large 
dots, squares, and open circles. 

used 128 3 zeroth- level grid in a computational box 
of 30/i _1 Mpc. The second set of simulations used 
256 3 grid in a 25ft. _1 Mpc box and had higher mass 
resolution. In the simulations used in this paper, 
the threshold for cell refinement (sec above) was 
low on the zeroth level: nthresh(O) — 2. Thus, ev- 



ery zeroth-level cell containing two or more par- 
ticles was refined. This was done to preserve 
all small-scale perturbations present in the ini- 
tial spectrum of perturbations. The threshold was 
higher on deeper levels of refinement. For the first 
set of simulations it was n t hrcsh = 2 at the first 
refinement level and n t hresh = 3 for all higher lev- 
els. For the second simulation the thresholds were 
"•thresh = 3 and n t hrosh = 4 for the first level and 
higher levels, respectively. 

In addition to effects introduced by limited 
mass and force resolution, integration errors of 
particle trajectories may affect the innermost re- 
gions of halos. The local dynamical time for par- 
ticles moving in these regions is quite short. For 
example, the period of a particle on a circular or- 
bit of radius l/i _1 kpc around the center of halo A 
is only 0.5% of the Hubble time. Therefore, if the 
time step is not sufficiently small, numerical errors 
in these regions will tend to grow. Even for small 
time steps errors exist and tend to alter the den- 
sity distribution in the centers of halos over some 
limited range of scales. 

All of our simulations were started at Zi = 60 
and the step in the expansion parameter was cho- 
sen to be Aao = 2 x I0~ 3 for particles located on 
the zeroth base grid. This gives about 500 steps 
for particles located in the zeroth level for an entire 
run to z = 0. We have done a test run with twice 
smaller time step for a halo of mass comparable 
(but with smaller number of particles) to the mass 
of halos studied in this paper. We did not find any 
significant differences in the resulting halo profile. 
For both sets of simulations, the highest level of 
refinement was ten for the largest mass resolution, 
which corresponds to 500 x 2 10 « 500, 000 time 
steps at the tenth refinement level. Some simula- 
tions were rerun with smaller number of particles. 
They did not reach the highest levels of refine- 
ment, and, thus, they had fewer steps. For exam- 
ple, halo D2 has reached only 7 levels of refinement 
and had only 500 x 2 7 « 64, 000 time steps. 

In the following sections we present density pro- 
files of four halos. The halo A was the only halo 
selected for re-simulation in the first set of simu- 
lations. It was relatively quiescent at z — and 
had no massive neighbors. The halo was located 
in a long filament bordering a large void and was 
about 10 Mpc away from the nearest cluster-size 
halo. After the high-resolution simulation was 



completed we found that the nearest galaxy-size 
halo was about 5 Mpc away. The halo had a fairly 
typical merging history with M(t) track slightly 
lower than the average mass growth predicted us- 
ing the extended Press-Schechter model. The last 
major merger event occurred at z ps 2.5; at lower 
redshifts the mass growth (the mass in this time 
interval has grown by a factor of three) was due 
to slow and steady mass accretion. 

The halos B, C, and D were identified in the 
second set of simulations and were selected among 
halos residing in a well defined filament. Two of 
the halos (B and C) are neighbors located about 
0.5 Mpc from each other. The third halo was 
2 Mpc away from this pair. Thus, the halos were 
not selected to be too isolated as was the case in 
the first set of runs. Moreover, the simulation was 
also analyzed at both z = and at z = 1 (when 
halos are more likely to be less relaxed). There- 
fore, the halo A can be considered as an example 
of a rather isolated well-relaxed halo. In many re- 
spects, this halo is similar to halos simulated by 
other research groups that used multiple mass res- 
olution techniques. The halos B, C, and D from 
the second set of simulations can be viewed as rep- 
resentative of more typical halo population located 
in more crowded environments. 

Parameters of the simulated dark matter halos 
are listed in Table 1. Columns in the table present 
(I) halo's "name" (halos Ai, A2, A3 are the halo A 
re-simulated three times with different mass and 
force resolutions); (2) redshift at which the halo 
was analyzed; (3) the number of particles within 
the virial radius; (4) the smallest particle mass 
in the simulation; (5) formal comoving force res- 
olution (cells size at the highest refinement level) 
achieved in the simulation. 

3. Results 

3.1. Comparison of the NFW and the 
Moore et al. profiles 

Before we fit analytical profiles to profiles of 
simulated dark matter halos or compare them to 
the observed rotation curves, it is instructive to 
compare different analytical approximations. Al- 
though the NFW and Moore et al. profiles predict 
different behavior of p(r) in the central regions of 
a halo, the scale at which this difference becomes 
significant depends on the specific values of the 
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Table 1 
Simulation parameters 



Halo 


2 




1 v part 


'"-part 


u tormal 


(1) 


(2) 




(3) 


(4) 


(5) 


Ai 





1 


.2 x 10 5 


1.6 x 10 7 


0.23 


A 2 





1 


.5 x 10 4 


1.3 x 10 8 


0.91 


A 3 





1 


.9 x 10 3 


1.1 x 10 9 


3.66 


Bi 





1 


.0 x 10 6 


1.2 x 10 6 


0.10 


B 2 





1 


.5 x 10 4 


1.2 x 10 6 


0.76 


B 3 


1 


7. 


.1 x 10 5 


1.2 x 10 6 


0.19 


Ci 





1 


.1 x 10 6 


1.2 x 10 6 


0.10 


c 2 





1 


.6 x 10 4 


7.7 x 10 7 


0.76 


c 3 


1 


5 


.0 x 10 5 


1.2 x 10 6 


0.19 


Di 





1 


.3 x 10 6 


1.2 x 10 6 


0.10 


D 2 





2 


.0 x 10 4 


7.7 x 10 7 


0.76 


D 3 


1 


7. 


.9 x 10 5 


1.2 x 10 6 


0.19 



halo's characteristic density and radius. Table 2 
presents the parameters and statistics associated 
with the two analytical profiles. For the NFW 
profile more information can be found in Klypin 
et al. (1998), Lokas & Mamon (2000), and Widrow 
(2000). 

Each profile is defined by two independent pa- 
rameters. We choose these to be the characteris- 
tic density p s and radius r s . In this case all ex- 
pressions describing the properties of the profiles 
have a simple form and do not depend on the con- 
centration. Both the concentration and the virial 
mass appear only in the normalization of the ex- 
pressions. The choice of the virial radius (e.g., 
Lokas & Mamon 2000) as a scale unit results in 
more complicated expressions with explicit depen- 
dence on the concentration. In this case, one has 
to be careful about the definition of the virial ra- 
dius, as there are several definitions in the liter- 
ature. For example, it is often defined as the ra- 
dius, r 20 o, within which the average density is 200 
times the critical density. In this paper the virial 
radius is defined as the radius within which the 
average density is equal to the density predicted 
by the top-hat model: it is 5th times the aver- 
age matter density in the Universe. In the case of 
fio = 0.3 models the virial radius defined in this 



way is about 30% larger than r 20 o (e.g., Eke et al. 
1998). 

There is no unique way of defining a consistent 
concentration for the different analytical profiles. 
Again, it is natural to use the characteristic radius 
r s to define the concentration: c = r v j r /r s . This 
simplifies the expressions. At the same time, if we 
fit the dark matter halo with the two profiles, we 
will get different concentrations because the values 
of the corresponding r s will be different. Alterna- 
tively, if we choose to match the outer regions of 
the profiles (say, r > r s ) as closely as possible, we 
may choose to change the ratio of the characteris- 
tic radii r S! NF\vAs, Moore hi such a way that both 
profiles reach the maximum circular velocity v c i IC 
at the same physical radius r max . In this case, the 
formal concentration of the Moore et al. profile is 
1.72 times smaller than that of the NFW profile. 
Indeed, with this normalization profiles look very 
similar in the outer parts as one finds in Figure 
3. Table 2 also gives two other "concentrations". 
The concentration C1/5 is defined as the ratio of 
virial radius to the radius, which encompasses 1/5 
of the virial mass (Avila- Reese et al. 1999). For 
halos with Cnfw ~ 5.5 this 1/5 mass concentra- 
tion is equal to Cnfw- One can also define the 
concentration as the ratio of the virial radius to 
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the radius at which the logarithmic slope of the 
density profile is equal to —2. This scale corre- 
sponds to r s for the NFW profile and w 0.35r s for 
the Moore et al. profile. 

Figure 3 presents the comparison between the 
analytic profiles normalized to have the same virial 
mass and the same radius r max . We show results 
for halos of low and high values of concentration 
representative of cluster- and low-mass galaxy ha- 
los, respectively. The bottom panels show the pro- 
files, while the top panels show the corresponding 
logarithmic slope as a function of radius. The fig- 
ure shows that the two profiles are very similar 
throughout the main body of the halos. Only in 
the very central region do the differences become 
significant. The difference is more apparent in the 
logarithmic slope than in the actual density pro- 
files. Moreover, for galaxy-mass halos the differ- 
ence sets in at a rather small radius < 0.01r v ; r , 
which would correspond to scales < 1 kpc for the 
typical dark matter dominated dwarf and LSB 
galaxies. At the observationally interesting scales 
the differences between NFW and Moore et al. 
profiles are fairly small and the NFW profile pro- 
vides an accurate description of the halo density 
distribution. 

Note also that for galaxy-size (e.g., high- 
concentration) halos the logarithmic slope of the 
NFW profile has not yet reached its asymptotic 
inner value of —1 even at scales as small as 
0.01r v ; r . At this distance the logarithmic slope 
of the NFW profile is ~ —1.4 — 1.5 for halos with 
mass ~ 10 12 ft.~ 1 M Q . For cluster-size halos this 
slope is « —1.2. This dependence of the slope at a 
given fraction of the virial radius on the virial mass 
of the halo is very similar to the results plotted 
in Figure 3 of Jing & Suto (2000). These authors 
interpreted it as evidence that halo profiles are 
not universal. It is obvious, however, that their 
results are consistent with NFW profiles and the 
dependence of the slope on mass can be simply a 
manifestation of the well-studied c v ; r (M) relation. 

The NFW and Moore et al. profiles can be 
compared in a different way. We can approximate 
the Moore et al. halo of a given concentration with 
the NFW profile. Fractional deviations of the fits 
depend on the halo concentration and on the range 
of radii used for the fits. A low-concentration halo 
has larger deviations, but even for C — 7 case, the 
deviations are less than 15% if we fit the halo at 
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Fig. 3. — Comparison of the Moore et al. and the 
NFW profiles. Each profile is normalized to have the 
same virial mass and the same radius of the maximum 
circular velocity. Left panels: High-concentration halo 
typical of small galaxy-size halos Cnfw = 17. Right 
panels: Low-concentration halo typical of cluster-size 
halos. The deviations are very small (< 3%) for radii 
r > r s /2. Top panels show the local logarithmic slope 
of the profiles. Note that for the high concentration 
halo the slope of the profile is significantly larger than 
the asymptotic value -1 even at very small radii r « 
0.01r vir . 

scales 0.01 < r/r v ; r < 1. For a high-concentration 
halo with C = 17, the deviations are much smaller: 
less than 8% for the same range of scales. 

To summarize, we find that the differences be- 
tween the NFW and the Moore et al. profiles are 
very small (Ap/p < 10%) for radii above 1% of 
the virial radius for typical galaxy-size halos with 
Cnfw ^ 12- The differences are larger for halos 
with smaller concentrations. In the case of the 
NFW profile, the asymptotic value of the central 
slope 7 = — 1 is not achieved even at radii as small 
as l%-2% of the virial radius. 

3.2. Convergence study 

The effects of numerical resolution can be stud- 
ied by resimulating the same objects with higher 
force and mass resolution and with a larger num- 
ber of time steps. In this study we performed 



8 



Table 2 

Comparison of NFW and Moore et al. profiles 



Parameter I NFW Moore et al. 



Density 

x = r/r s 


n- P g 
P x{l + xf 

p oc x~ 3 for x ^ 1 
p oc x -1 for x <C 1 
p/ps = 1/4 at x = 1 


P x^l+x) 1 - 5 

p oc x -3 for x ^> 1 
p oc x -1 ' 5 for x <C 1 
p/p s = 1/2 at x = 1 


Mass 

M = 4irp s r 3 s f(x) 
= M vi J(x)/f(C) 

M vil = ^Pcr^O^top-hatrvir 


/(x)=ln(l + x)- T ^ 


/(x) = fln(l + x 3 / 2 ) 


Concentration 
C = r vil /r s 


CnFW = l-72CMooro 

for halos with the same M v i T and r max 
/-i . ~ Cnfw 


CMoorc = CnFw/1-72 
r< _ Cmooto 


V5 ~ o.86/(C NF w) +0.1363 
error less than 3% for Cnfw =5-30 

C 7= _2 = Cnfw 


1/5 id +Cj 1/5 - 1] 2/3 

^ CMoorc 
~ lc 3/W _ -,2/3 
[°Moorc X J 

C 7= -2 = 2 3 ' 2 CMooro 

~ 2.83CMooro 


Circular Velocity 

2 GM vir C f(x) 

Vcilc r vir x/(C) 

^2 *^max f\pt) 

2 GM vir 

v ■ = 

™ „ . 

1 vir 


«£max ~ 2.15 

c 

<ax~0.216^. r ^. ( ^ 
p/p s w 1/21.3 at x = 2.15 


^max ^ 1-25 

v 2 max « 0.466^ 

p/p s w 1/3.35 at x = 1.25 
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simulations of the same halos with increasingly 
higher mass resolution. In the ART code simu- 
lations the subsequent mesh refinements are done 
when particle density in a mesh cell exceeds a spec- 
ified threshold. The mass resolution is thus tightly 
linked with the highest achievable spatial resolu- 
tion. 

Table 3 gives parameters of the halos and pa- 
rameters of their fits. The first and the second 
columns give the halo name (Table 1) and the red- 
shift at which the halo was studied. Columns (3-5) 
present virial mass, radius, and the maximum cir- 
cular velocity of the halo. Columns (6-8) present 
parameters of the fits: the halo concentration as 
estimated using the NFW profile and the maxi- 
mum relative errors of the NFW and Moore et al 
fits. The bottom panel in figure 4 shows density 
profiles for the simulations of halo A (see Table 1). 
Here, as in the Fig. la in Moore et al. (1998), all 
profiles are plotted down to the formal force reso- 
lution of the corresponding run. Although it may 
appear that density profiles have not converged in 
the central region and that the low resolution sim- 
ulations produce erroneous results, this is simply 
an artifact of plotting the profiles below the actual 
numerical resolution (or convergence scale). The 
top panel in figure 4 and figure 5 show profiles of 
halos A, B, C, and D plotted down to four formal 
resolutions of the simulations (4 mesh cells at the 
highest refinement level). The figures show that in 
this case density profiles in lower resolution sim- 
ulations arc in very good agreement with profiles 
in high-resolution runs at all radii. For example, 
there are no systematic differences in the logarith- 
mic slope of the profile at a given distance and 
we find no significant change in the concentration 
parameters or the maximum circular velocities of 
halos (see Table 3). 

Closer examination of the bottom panel in Fig- 
ure 4, shows that profiles have not converged at 
two formal resolutions. This is at odds with our 
convergence study in Kravtsov et al. (1998). We 
attribute this difference to the fact that mass reso- 
lution in the latter study was kept fixed when force 
resolution was varied. Our highest resolution run 
Ai , if considered including scales larger than two 
formal resolutions, is consistent with conclusion 
about the shallow central slope made in Kravtsov 
et al. (1998). Indeed, if the profiles are consid- 
ered down to the scale of two formal resolutions (a 
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Fig. 4. — Density profiles of halo A simulated with 
different mass and force resolutions. Bottom panel: 
The profiles are plotted down to the formal force res- 
olution of each simulation. Because the plot shows 
results below the actual force resolution, one gets a 
wrong impression that the profile gets steeper and the 
concentration increases when the mass resolution is in- 
creased. Top panel: the profiles plotted down to four 
formal resolutions. For vastly different mass and force 
resolutions the convergence is reached at these scales. 
At the scale of a few percent of the virial radius the 
density profile is visibly steeper than the limiting slope 
7 = — 1 of the NFW profile. This is consistent with 
the NFW profile for a halo of this concentration. 

scale smaller than the smallest converged scale), 
the density profile slope in the very central part 
of the profile r < 0.01r v i r is close to 7 = —0.5. In 
light of the results shown in Figure 4, it is clear 
that this is an artifact of underestimating true con- 
vergence scale and conclusions about shallow cen- 
tral density distributions made in Kravtsov et al. 
(1998) are therefore incorrect. 
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Table 3 
Halo parameters 



Halo 


2 


M vir 


D 




I i "NT T? \ S T 

^ IN r W 


RelErr 


RelErr 






h~ 




h~ "^kpc 


km J s 




NFW 




(1) 


(2) 




(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


Ai 





2.0 


x 10 12 


257 


247.0 


17.4 


0.17 


0.20 


A 2 





2.1 


x 10 12 


261 


248.5 


16.0 


0.13 


0.16 


A 3 





2.0 


x 10 12 


256 


250.5 


16.6 


0.16 


0.10 


Bi 





1.2 


x 10 12 


215 


199.5 


15.6 


0.30 


0.14 


B 2 





1.1 


x 10 12 


213 


205.0 


16.5 


0.15 


0.14 


B 3 


1 


8.5 


x 10 11 


241 


195.4 


12.3 


0.23 


0.16 


Ci 





1.3 


x 10 12 


225 


190.6 


11.2 


0.29 


0.23 


c 2 





1.2 


x 10 12 


220 


184.9 


9.8 


0.11 


0.12 


c 3 


1 


6.8 


x 10 11 


208 


165.7 


11.9 


0.37 


0.20 


Di 





1.5 


x 10 12 


235 


213.9 


11.9 


0.15 


0.68 


D 2 





1.5 


x 10 12 


234 


216.8 


13.4 


0.10 


0.09 


D 3 


1 


9.6 


x 10 11 


245 


202.4 


9.5 


0.25 


0.60 



It is not clear which numerical effect determine 
the convergence scale. It is likely that this scale is 
determined by a complex interplay of all numerical 
effects. For the ART simulations we found empir- 
ically that the scale above which the density does 
not deviate (deviations were less than 10%) from 
results of higher resolution simulation is > 4 for- 
mal force resolutions or containing more than 200 
particles, whichever is larger. The limit very likely 
depends on particular code used and is not univer- 
sal. Figure 4, top panel, shows that for the halo A, 
convergence for vastly different mass and force res- 
olution is reached for scales > 4 formal force reso- 
lutions (all profiles in this figure are plotted down 
to the radius of 4 formal force resolutions). For 
all resolutions, there are more than 200 particles 
within the radius of four resolutions from the halo 
center. For the highest resolution simulation (halo 
Ai) convergence is reached at scales > 0.005r v i r , 
assuming convergence at 4 times the formal reso- 
lution as found for halos A 2 and A 3 . For halos B, 
C, D (figure 5) this criterion also worked, but was 
mostly defined by the number of particles (more 
than 200-300 particles for convergence). 



3.3. Halo profiles 

In order to judge which analytical profile pro- 
vides a better description of the simulated profiles 
we fitted the NFW and Moore et al. analytic pro- 
files. Figure 6 presents results of the fits for halo 
A and shows that both profiles fit the simulated 
profile equally well: fractional deviations of the 
fitted profiles from the numerical one are smaller 
than 20% over almost three decades in radius. It is 
thus clear that the fact that the numerical profile 
has slope steeper than —1 at the scale of ~ 0.01r v ir 
does not mean that a good fit of the NFW profile 
(or even analytic profiles with shallower asymp- 
totic slopes) cannot be obtained. Figure 7 shows 
the fitting of halos in the second set of simulations. 
Each halo in the plot has more than a million par- 
ticles - ten times more than halo A. One would 
naively expect that this increase in the resolution 
should clearly show which profile makes a better 
fit. Indeed, more particles and better resolution 
gave smaller deviations, but the fits became bet- 
ter for both approximations. For example, at 1% of 
the virial radius of the halo D the deviations were 
3.6% for the NFW profile and 6.2% for the Moore 
et al. profile - down from 20% for the halo A at the 
same distance. The Moore et al. approximation 
gave a better fit for halos B and C, but not for halo 



11 




1 5 10 50 100 

R(kpc/h) 



Fig. 5. — Convergence of profiles for halos B, C, D. 
All the curves in this plot start at radii, which con- 
tain at least 200 particles and which are larger than 
4 formal resolutions. For clarity, the profiles of ha- 
los C and B were shifted downwards by factors of 10 
and 100, respectively. The full curves are for halos 
Bi, Ci, Di, which have more than a million particles 
each. The dashed curves show profiles of halos B2, 
C 2 , D 2 with (15 - 20) x 10 3 particles. The error bars 
indicate 2<r fluctuations of density at each radius due 
to moving satellites and residual oscillations inside the 
halos. The dispersion was estimated using density pro- 
files at four time moments in the interval 2 = — 0.03 
(Ri4x 10 8 yrs). 

D. The NFW approximation was less accurate on 
intermediate scales around (0.03 — 0.1)i? v i r , but 
the errors were quite small. Thus, both approxi- 
mations gave comparable results. 

There is definitely a certain degree of degener- 
acy in fitting various analytic profiles to numerical 
results. Figure 8 illustrates this further by show- 
ing results of fitting profiles (solid lines) of the 
form p(r) oc (r/r )^[l + (r /r )"]- (/3 - Q)/7 to the 
same simulated halo profile (halo Ai) shown as 
solid circles. The legend in each panel indicates 
the corresponding values of a, (3, and 7 of the fit; 
the digit in parenthesis indicates whether the pa- 
rameter was kept fixed (0) or not (1) during the 
fit. The two right panels show fits of the NFW and 
Moore et al. profile; the bottom left panel shows 
fit of the profiles used by Jing & Suto (2000). The 
top left panel shows a fit in which the inner slope 
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Fig. 6. — Fits of the NFW and Moore et al. halo 
profiles to the profile of halo Ai (bottom panel). The 
top panel shows fractional deviations of the analytic 
fits from the numerical profile. Note that both ana- 
lytic profiles fit numerical profile equally well: frac- 
tional deviations are smaller than 20% over almost 
three decades in radius. 

was fixed but a and (3 were fit. The figure shows 
that all four analytic profiles can provide a good 
fit to the numerical profile in the whole range of 
resolved scales: 0.005 — lr v ; r . 

As we mentioned in § 2.3, the halo A analyzed 
in the previous section is somewhat special be- 
cause it was selected as an isolated relaxed halo. 
Halos, which are not very isolated and relaxed are 
also interesting. After all, they represent the ma- 
jority of all halos. When we compare observed 
rotation curves of galaxies with predicted circular 
velocity curves of the halos, we do not know if the 
galaxy host halo is well relaxed or not. In order to 
reach unbiased conclusions, we will present anal- 
ysis of halos from the second set of simulations 
at redshifts z = and z = 1 (halos Bi^ 3 , Ci,3, 
and D13), which were not selected to be relaxed 
or isolated. Note that these halos did not have 
major mergers immediately prior to the epoch of 
analysis, which could produce large distortions of 



12 



«(r)«(r/r )-*[l+(r/r )«]-«<-°>/7 



10" 
10 s 

K3 10< 

+ 



' 1 1 1 

|~»-^ halo 1 

• 1.0 (0) -j 
I a- 0.55 (1) J 

: /s= 4.03 (i) ^1 

1" r G = 0.23 j 


:^«^ halo l : 

r *> i 

r 7= 1.0 (0) - 
: _ a= 1.0 (0) J 

(S= 3.o (o) 

r r, - 0.06 (NFW) 


•^H h] l-H 1 hHffl 

1 halo 1 1 

f 7= 1-5 (0) \^ 

L a= 1.0 (0) J 
: fi 3.0 (0) V ! 

f r = 0.120 (ling & Suto) I 




^1 1 "Ml 1 1 1 1 '1' 1 1 1 1 III 

1 halo 1 | 

r y i-5 (o) >i 

1 <* = 1.5 (0) J 
E P- 3.0 (0) ^E 

f r = 0.120 (Moore) 





0.01 0.1 1 0.01 0.1 



r/r . r/r . 

Fig. 8. — Analytic fits to the density profile of the 
halo Ai from our set of simulations. The fits are of 
the form p(r) oc (r/r )~ 7 [l + (r /r ) a ]- {f) - a)h . The 
legend in each panel indicates the corresponding val- 
ues of a, (3, and 7 of the fit; the digit in parenthesis 
indicates whether the parameter was kept fixed (0) or 
not (1) during the fit. Note that various sets of param- 
eters a, (3, 7 provide equally good fits to the simulated 
halo profile in the whole range resolved range of scales 
« 0.005 — lr v ir- This indicates a large degree of de- 
generacy in parameters a, [3, and 7 

their profiles. Based on the results of the conver- 
gence study presented in the previous section, we 
will consider profiles of these halos only at scales 
above four formal resolutions and not less than 200 
particles. There is an advantage in analyzing halos 
at a relatively high redshift. Halos of a given mass 
will have lower concentration (see Bullock et al. 
2000). Lower concentration implies a large scale 
at which the asymptotic inner slope is reached. 

We found that substantial substructure is 
present inside the virial radius in all three ha- 
los at z = 1. Figure 7 shows profiles of these 
halos at z = (top) and z = 1 (bottom). The 
z = profiles are smoother than profiles at z = 1. 
Note that bumps and depressions visible in the 
profiles have amplitude that is significantly larger 
than the shot noise. Halo C3 appeared to be the 
most relaxed of the three halos. This halo had 
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Fig. 9. — Circular velocity profiles for the halos Bi, 
Ci, and Di normalized to halo's virial velocity. Halos 
are well resolved on all shown scales. Although the 
halos have very similar masses, the profiles are very 
different; the differences are due to real differences in 
the concentration parameters. 

its last major merger somewhat earlier than the 
other two. Halo D3 had a major merger event at 
z w 2. A remnant of the merger is still visible as a 
bump at r ~ 100/i _1 kpc. The non-uniformities of 
profiles caused by substructure may substantially 
bias analytic fits if one uses the entire range of 
scales below the virial radius. Therefore, we used 
only the central, presumably more relaxed, regions 
in the analytic fits: r < 50/i _1 kpc for halo D and 
r < 100/i -1 kpc for halos B and C (fits using only 
central 50ft. _1 kpc did not change results). 

The best fit parameters were obtained by min- 
imizing the maximum fractional deviation of the 
fit: max[abs(log ( Ofit — logphalo)]- Minimizing the 
sum of squares of deviations (x 2 ), as is often done, 
can result in larger errors at small radii with the 
false impression that the fit fails because it has 
a wrong central slope. The fit that minimizes 
maximum deviations improves the NFW fit for 
points in the range of radii (5 — 20)/i~ 1 kpc, where 
the NFW fit would appear to be below the data 
points if the fit was done by the x 2 minimization. 
For example, if we fit halo B by minimizing x 2 > 
the concentration slightly decreases from 12.3 (see 
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Fig. 7. — Fits of the NFW (dotted curves) and Moore et al. (dot-dashed curves) profiles to the density distributions 
of halos B, C, D (solid curves) in the highest resolution run at z = (top left). For clarity, the profiles of halos C 
and B were shifted down by factors of 10 and 100, respectively. Error bars present 2a fluctuations of density profiles 
at different time moments estimated over the last ~ 4 x 10 s yr of the run. Top right panel shows fractional deviations 
of the fitted profile from the simulated halo profile as a function of scale. The bottom panels show the same but for 
the halos B, C, and D at z = 1. Note that halos at this redshift are not yet well relaxed as indicated by the density 
fluctuations in the peripheral parts of the profiles. 

Table 1) to 11.8, the maximum error slightly in- points for most of the points at small radii, 

creases to 27%, but the fit goes below the data We have also fitted density distribution of halo 
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B assuming even more stringent limits on the ef- 
fects of numerical resolution. We fitted the halo 
starting at the scale equal to six times the for- 
mal resolution, minimizing the maximum devia- 
tion. Inside this radius there were about 900 par- 
ticles. Resulting parameters of the fit were close 
to those in Table 1: Cnfw = 11.8, and maximum 
error of the NFW fit was 17%. 

We found that for halos B and C the errors in 
the Moore et al. fits were systematically smaller 
than those of the NFW fits, though the differences 
were not dramatic. But Moore ct al. fit poorly in 
the case of halo D. It formally gave very small 
errors, but at the expense of unreasonably small 
concentration Cnfw = 2. When we constrained 
the approximation to have about twice larger con- 
centration as compared with the best NFW fit, 
we were able to obtain a reasonable fit (this fit 
is shown in Figure 7). Nevertheless, the central 
density distribution is fit poorly in this case. 

Therefore, our analysis does not show that one 
analytic profile is better then the other for descrip- 
tion of the density distribution in simulated halos. 
Despite the larger number of particles per halo and 
the lower concentrations of z = 1 halos, results are 
still inconclusive. The Moore at al. profile is a bet- 
ter fit to the profile of halo C; the NFW profile is 
a better fit to the central part of the halo D. Halo 
B represents an intermediate case where both pro- 
files provide equally good fits (similar to the anal- 
ysis of halo A). Remarkably, the same conclusions 
hold for the halo profiles at z = 0. 

Both at z = and z — 1, there are real de- 
viations in parameters of halos of the same mass. 
We find the same differences in estimates of Ci/ 5 
concentrations, which do not depend on specifics 
of an analytic fit. The central slope at around 
1 kpc also changes from halo to halo. Halos B and 
C have the same virial radii and nearly the same 
circular velocities, yet their concentrations are dif- 
ferent by 30%. Indeed, halos in the Table 3 have 
similar masses in the range (1.2 — 2) x 10 12 /i _1 M Q . 
If halos had a universal profile - a shape, which 
depends only on halo mass, then we should expect 
that the circular velocity curves are very similar 
for our halos. Figure 9 shows circular velocities 
for halos B, C, and D, which have only 25% de- 
viations in their virial mass. The halos clearly do 
not have a universal one-parameter shape. There 
are substantial variations in the curves, which oc- 



cur at relatively large radii (~ 0.1 — 0.3i? v ; r ). The 
variations are due to differences in halo concentra- 
tion - each curve is well described by a NFW or 
Moore et al. profile but their concentrations are 
somewhat different. Our three halos clearly con- 
stitute a small sample. Bullock et al. (2000) and 
Jing (2000) studied the spread of halo concentra- 
tions in a large sample of halos. For a given mass 
it was found that halos have 20-50% variations in 
the concentration at la level, which is consistent 
with what we find for our halos. 

4. Discussion and conclusions 

We have analyzed a series of simulations with 
vastly different mass and force resolutions with the 
goal of studying density distribution in the cen- 
tral regions of galaxy-size dark matter halos. We 
used multiple mass simulations performed using 
the ART code; the simulations were performed 
with variable mass, force and temporal resolu- 
tions. In the highest resolution runs, we achieved a 
(formal) spatial dynamical range of 2 18 = 262144; 
the simulation was run with 500,000 steps for par- 
ticles at the highest level of refinement. 

Using these simulations, we have studied con- 
vergence of halo density profiles for different mass 
and force resolutions. We show that the halo pro- 
files converge at scales larger than a certain (true 
numerical resolution) scale defined by numerical 
effects. This scale is probably code dependent, but 
can be found for any numerical code by a conver- 
gence study. For the ART simulations presented 
here, the density profiles converged at the scale 
of four times the formal force resolution or the ra- 
dius containing more than 200 particles, whichever 
is larger. In this sense, our results are consistent 
with results of the "Santa Barbara" cluster com- 
parison project (Frenk et al. 1999): the density 
profiles of a cluster-size halo simulated with dif- 
ferent numerical codes and resolutions agree with 
each other at all resolved scales. 

In KKBP98 we have discussed convergence 
tests in which we varied force resolution while 
keeping mass resolution fixed. Using these tests we 
concluded that density profiles converge at radii 
twice the local formal resolution of a simulation. 
The convergence tests presented in this paper, 
however, show that mass resolution places more 
stringent conditions on the trustworthy range of 
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scales. Although we can reproduce our previous 
results (shallower than —1 density profiles at radii 
of two formal resolution), our new convergence 
tests show that these results were affected by lim- 
ited mass resolution. We conclude that we overes- 
timated our force resolution KKBP98 and that the 
conclusions about the shallow central slopes pre- 
sented there were an artifact of that overestimate. 
Other results in KKBP98 that focus on general 
halo characteristics, such as the scatter in profile 
shapes, and the agreement between the V max - 
Tmax relations of simulated dark halos and dark 
matter dominated dwarf and LSB galaxies, are 
valid. 

At scales above four times the formal resolution 
and containing more than 200 particles results pre- 
sented in this paper agree well with previous sim- 
ulations. For example, the concentration param- 
eter for the halos are in good agreement with the 
concentration-mass dependence (c(M)) presented 
in Bullock et al. (2000) based on the previous 
simulations. 

We can also reproduce results of convergence 
studies by Moore et al. (1998) and Ghigna et 
al. (1999). At first glance it may appear that 
our conclusions are in direct contradiction with 
these studies that concluded that at least sev- 
eral millions of particles arc needed to resolve a 
halo profile properly. However, the contradiction 
is only in interpretation rather than in the results 
themselves. For example, Figure 2 in Ghigna et 
al. shows a cluster profile simulated at 3 different 
mass and force resolutions. The conclusion the au- 
thors make based on this convergence study is (at 
least qualitatively) similar to our conclusion: the 
profiles converge at all mass resolutions at scales 
above 6e, where e is the spline force softening of 
their code 4 (roughly equivalent to our formal res- 
olution) . Why is this criterion is more severe than 
ours? The obvious reasons are differences in the 
force shape and other differences between numer- 
ical codes. However, it appears that the soften- 
ing in these simulations was set too low (the pro- 
files were "overresolved" ) and convergence scale 
is determined by the number of particles crite- 
rion rather than by force resolution. Indeed, we 

4 Note that profiles of lower resolution simulations perfectly 
agree with profile of the highest resolution run at scales 
above 3 — 4e. 



can estimate the number of particles within the 
softening scale in these simulations as N p (< e) w 
(47r/3)(/9 cr it/m p )(l + <5)e 3 , where e is the soften- 
ing length, m p is the particle mass, p cr n is the 
critical density, and S is the overdensity reached 
at the softening scales. For the HIRES simula- 
tion in Ghigna et al., e = 0.5ft. -1 kpc and m p = 
5.37 xlO 7 /^ 1 Mq and J ~ (0.5-=- 1.0) x 10 7 , which 
gives N p (< e) ~ 10 - 30. For their LOWRES sim- 
ulation this number is even lower ~ 5 — 10. The 
scale containing > 200 — 300 particles should be 
close to the convergence that is found in this study 
(e.g., 3 - 4e for the LOWRES run). 

In light of these considerations, the profiles in 
Fig 2 of Moore et al. (1998) are perfectly consis- 
tent with each other if considered at scales > 4e. 
Their results are thus perfectly consistent with our 
results and conclusions. One needs more parti- 
cles only if the softening is set too small and the 
inner regions are over-resolved. In this case the 
true resolution is set by the radius that contains at 
least a couple hundred particles. The main point 
is that the higher mass resolution is needed to 
probe deeper into the inner regions of halos. How- 
ever, if one is interested only in the profiles at, say, 
r > 0.02R vir (sufficient to determine halo's con- 
centration, maximum circular velocity, etc.) than, 
as Figure 2 in Ghigna et al. (1999) and Figure 5 
in this paper clearly show, ~ (2 — 5) x 10 5 within 
the virial radius is adequate. 

In this paper, we present results for halos that 
contain 1.2 x 10 5 to 1.6 x 10 6 particles within their 
virial radius. We show that for the galaxy-size ha- 
los (M vir = 7 x 10 11 /i- 1 M Q - 2 x lO^fr^M© and 
C = 9-17) both the NFW profile p oc r" 1 (l+r)~ 2 
and the Moore et al. profile p oc r _1 ' 5 (l + r 1 ' 5 ) -1 
provide fairly good fits of the simulated profiles 
with deviations of about 10% for radii larger than 
1% of the virial radius. For dwarf and LSB galax- 
ies commonly used for comparisons with model 
predictions, this corresponds to scales lkpc. 
Therefore, the debate about which analytic pro- 
file provides a better description of the CDM halo 
profiles may well be irrelevant for comparisons to 
measured galaxy rotation curves. Such compar- 
isons are also subject to other uncertainties, one 
of which is limited spatial extent of the observed 
rotation curves. The particular shape of the inner 
density distribution may be important for galaxy 
cluster observations, however. Cluster-size ha- 
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los are predicted to have smaller concentrations, 
which means that the scale at which differences be- 
tween the NFW and Moore et al. profiles become 
significant is larger and is observationally relevant. 

These results are consistent with results of 
Moore et al. (1999) who found that galaxy-size 
halos in CDM models are well described by the 
p oc r _1 ' 5 (l + r 1,5 ) -1 profile. The authors used 
simulations of the standard CDM model with mass 
and force resolution similar to our simulations 
(2 x 10 6 M Q and 1 kpc, respectively). They found 
that the NFW profile that was fitted at radii above 
3% of the halo's virial radius, underpredicts the 
density at smaller radii by up to 20 — 30%. This is 
consistent with out results for halos B and C, al- 
though not for halo D, which probably indicates a 
certain degree of variance among profiles of halos 
of the same mass (Jing 2000; Bullock et al. 2000; 
Avila-Reese et al. 1999). 

Jing & Suto (2000) simulated formation of 
galaxy-, group-, and cluster-size halos with res- 
olution similar to our highest-resolution simula- 
tions. They also found that the NFW profile fit 
to the outer regions (radii above few percent of 
the virial radius) underestimates the density in 
the innermost regions of the halo. The degree of 
discrepancy, however, appears to be different for 
different halos (sec their Fig. 2), which is consis- 
tent with our conclusions. The authors conclude 
that the shapes of the density profiles vary from 
galaxy- to group-, to cluster-size halos. However, 
as we argued in § 3.1, their results could be inter- 
preted as the manifestation of concentration-mass 
relation instead. 

Most recently, Fukushige & Makino (2001) sim- 
ulated 12 halos with masses ranging from 4.3 x 
10 11 Mq to 7.9 x 10 14 M Q , again with resolution 
similar to the resolution of simulations presented 
here. The dependence of the inner logarithmic 
slope on the halo mass observed by Jing & Suto 
(2000) was not found in these simulations. In- 
stead, the steepest slope of the density profiles 
was found to be close to —1.5 for all masses. 
This is surprising because such dependence is ex- 
pected from the concentration-mass correlation 
(e.g., NFW; Bullock et al. 2000). It is not clear 
what could explain this discrepancy. Note also 
that setup of these simulations was somewhat dif- 
ferent from the setup that is usually used: the 
simulations followed halo collapse from a spheri- 



cally symmetric configuration centered on a gaus- 
sian density peak with no external tidal field in- 
cluded. 

We show that density profiles of halos that are 
not fully relaxed may contain real non-uniformities 
due to substructure and differences in the density 
distributions. These non-uniformities affect the fit 
quality for a particular analytic profile and result 
in somewhat different values of fitted parameters 
(e.g., concentration). This was clearly seen at red- 
shift z = 1 for halos B,C,D (Figure 7), which by 
that time have not yet relaxed. At redshift z = 
the halos were much more quiet and the deviations 
from the fits were much smaller. It is interesting to 
note that the non-equilibrium effects do not qual- 
itatively change the shape of the central parts of 
density profiles. For example, we found that at 
both redshifts the profile of halo C is best fit by 
the Moore et al. profile. The density profile of halo 
D, however, is best fit by the NFW profile again at 
both redshifts. Note that these three halos were 
simulated with the same mass and force resolution 
(indeed, in the same simulation). It seems that the 
main difference between these halos is their merger 
histories. We conclude, therefore, that differences 
in merger history and/or different degree of sub- 
structure in halos of the same mass may explain 
the scatter in profile shapes and concentration pa- 
rameters found in previous studies (KKBP98; Jing 
2000; Bullock et al. 2000). 

In view of the current less confusing situation 
regarding the theoretical predictions for CDM pro- 
files, it is interesting to discuss how the theory 
compares to observations. Rotation curves of a 
number of dwarf and LSB galaxies have recently 
been re-examined using Ha observations and/or 
including corrections for beam-smearing in HI ob- 
servations (e.g., Swaters, Madore & Trcwhclla 
2000; van den Bosch et al. 2000). The results show 
that for the majority of galaxies, the Ha rotation 
curves are significantly different in their central 
regions than the rotation curves derived from HI 
observations. This indicates that the HI rotation 
curves are affected by beam smearing (Swaters et 
al. 2000). This also implies that beam smearing 
may be at least partly responsible for the univer- 
sal shape of the LSB rotation curves discussed in 
KKBP98. 

It is possible that part of the discrepancy be- 
tween the rotation curves of the different trac- 
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ers may be due to real differences in the kine- 
matics of the two gas components (ionized and 
neutral hydrogen). Preliminary comparisons be- 
tween the new Ha rotation curves and model pre- 
dictions show that NFW density profiles arc con- 
sistent with the observed shapes of the rotation 
curves (van den Bosch & Swatcrs 2000; Navarro 
& Swaters 2000). Moreover, cuspy density profiles 
with inner logarithmic slopes as steep as ~ —1.5 
also seem to be consistent with the data (van den 
Bosch & Swaters 2000). A separate concern is not 
the shape of the inner density profile, but rather 
the value of the central density. There are in- 
dications that CDM halos are too concentrated 
(Navarro & Swaters 2000; McGaugh et al. 2000; 
Navarro & Stcinmetz 2000; Firmani et al. 2000) 
in comparison with galactic halos. However, van 
den Bosch & Swaters (2000) have argued, based 
on detailed modeling of adiabatic contraction and 
beam-smearing, that dwarf galaxy concentrations 
are in fact consistent with the observed distribu- 
tion in ACDM halos. Thus, although the shape 
of galactic rotation curves may be not as differ- 
ent from predictions as was thought before, the 
halo concentrations derived from observations are 
alarmingly low. The recent observational progress 
made in this field promises to resolve and clar- 
ify this issue in the near future. At larger scales 
(r > 100 kpc), the constraints from weak galaxy- 
galaxy lensing (Fischer et al. 2000; Smith et al. 
2000) should be very useful in constraining the 
overall profile and concentrations of galactic ha- 
los. 

In summary, the study presented here is aimed 
to clarify the issue of convergence of the density 
profiles of CDM halos. We show that convergence 
can be reached regardless of the mass resolution, 
although the convergence scale does depend on the 
mass resolution: the higher mass resolution results 
in smaller convergence scale for the same objects, 
but it does not affect the outer parts of the profile. 
Our results also indicate that there is a real scatter 
in shapes of density profiles and halo parameters. 
Larger systematic studies currently underway may 
put these conclusions on a firmer footing. 
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